In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
from construct import *
In [2]:
frame_size = 223*5
frames = np.fromfile('bepi_colombo_frames.u8', dtype = 'uint8')
frames = frames[:frames.size//frame_size*frame_size].reshape((-1, frame_size))
frames.shape[0]
Out[2]:
CCSDS TM frame CRC is CRC16_CCITT_FALSE from this calculator.
In [3]:
crc_table = [0x0000, 0x1021, 0x2042, 0x3063, 0x4084, 0x50a5, 0x60c6, 0x70e7,
0x8108, 0x9129, 0xa14a, 0xb16b, 0xc18c, 0xd1ad, 0xe1ce, 0xf1ef,
0x1231, 0x0210, 0x3273, 0x2252, 0x52b5, 0x4294, 0x72f7, 0x62d6,
0x9339, 0x8318, 0xb37b, 0xa35a, 0xd3bd, 0xc39c, 0xf3ff, 0xe3de,
0x2462, 0x3443, 0x0420, 0x1401, 0x64e6, 0x74c7, 0x44a4, 0x5485,
0xa56a, 0xb54b, 0x8528, 0x9509, 0xe5ee, 0xf5cf, 0xc5ac, 0xd58d,
0x3653, 0x2672, 0x1611, 0x0630, 0x76d7, 0x66f6, 0x5695, 0x46b4,
0xb75b, 0xa77a, 0x9719, 0x8738, 0xf7df, 0xe7fe, 0xd79d, 0xc7bc,
0x48c4, 0x58e5, 0x6886, 0x78a7, 0x0840, 0x1861, 0x2802, 0x3823,
0xc9cc, 0xd9ed, 0xe98e, 0xf9af, 0x8948, 0x9969, 0xa90a, 0xb92b,
0x5af5, 0x4ad4, 0x7ab7, 0x6a96, 0x1a71, 0x0a50, 0x3a33, 0x2a12,
0xdbfd, 0xcbdc, 0xfbbf, 0xeb9e, 0x9b79, 0x8b58, 0xbb3b, 0xab1a,
0x6ca6, 0x7c87, 0x4ce4, 0x5cc5, 0x2c22, 0x3c03, 0x0c60, 0x1c41,
0xedae, 0xfd8f, 0xcdec, 0xddcd, 0xad2a, 0xbd0b, 0x8d68, 0x9d49,
0x7e97, 0x6eb6, 0x5ed5, 0x4ef4, 0x3e13, 0x2e32, 0x1e51, 0x0e70,
0xff9f, 0xefbe, 0xdfdd, 0xcffc, 0xbf1b, 0xaf3a, 0x9f59, 0x8f78,
0x9188, 0x81a9, 0xb1ca, 0xa1eb, 0xd10c, 0xc12d, 0xf14e, 0xe16f,
0x1080, 0x00a1, 0x30c2, 0x20e3, 0x5004, 0x4025, 0x7046, 0x6067,
0x83b9, 0x9398, 0xa3fb, 0xb3da, 0xc33d, 0xd31c, 0xe37f, 0xf35e,
0x02b1, 0x1290, 0x22f3, 0x32d2, 0x4235, 0x5214, 0x6277, 0x7256,
0xb5ea, 0xa5cb, 0x95a8, 0x8589, 0xf56e, 0xe54f, 0xd52c, 0xc50d,
0x34e2, 0x24c3, 0x14a0, 0x0481, 0x7466, 0x6447, 0x5424, 0x4405,
0xa7db, 0xb7fa, 0x8799, 0x97b8, 0xe75f, 0xf77e, 0xc71d, 0xd73c,
0x26d3, 0x36f2, 0x0691, 0x16b0, 0x6657, 0x7676, 0x4615, 0x5634,
0xd94c, 0xc96d, 0xf90e, 0xe92f, 0x99c8, 0x89e9, 0xb98a, 0xa9ab,
0x5844, 0x4865, 0x7806, 0x6827, 0x18c0, 0x08e1, 0x3882, 0x28a3,
0xcb7d, 0xdb5c, 0xeb3f, 0xfb1e, 0x8bf9, 0x9bd8, 0xabbb, 0xbb9a,
0x4a75, 0x5a54, 0x6a37, 0x7a16, 0x0af1, 0x1ad0, 0x2ab3, 0x3a92,
0xfd2e, 0xed0f, 0xdd6c, 0xcd4d, 0xbdaa, 0xad8b, 0x9de8, 0x8dc9,
0x7c26, 0x6c07, 0x5c64, 0x4c45, 0x3ca2, 0x2c83, 0x1ce0, 0x0cc1,
0xef1f, 0xff3e, 0xcf5d, 0xdf7c, 0xaf9b, 0xbfba, 0x8fd9, 0x9ff8,
0x6e17, 0x7e36, 0x4e55, 0x5e74, 0x2e93, 0x3eb2, 0x0ed1, 0x1ef0]
def crc16_ccitt_false(data):
crc = 0xffff
for d in data:
tbl_idx = ((crc >> 8) ^ d) & 0xff
crc = (crc_table[tbl_idx] ^ (crc << 8)) & 0xffff
return crc & 0xffff
In [4]:
crc_ok = np.array([crc16_ccitt_false(f) for f in frames]) == 0
In [5]:
np.average(crc_ok)
Out[5]:
In [6]:
np.sum(~crc_ok)
Out[6]:
In [7]:
plt.figure(figsize = (15,15), facecolor = 'w')
plt.imshow(frames[crc_ok,:], interpolation = 'nearest', aspect = 0.4)
Out[7]:
In [8]:
plt.figure(figsize = (30,15), facecolor = 'w')
plt.imshow(frames[crc_ok,:16], interpolation = 'nearest', aspect = 0.1)
Out[8]:
In [9]:
plt.figure(figsize = (30,15), facecolor = 'w')
plt.imshow(frames[crc_ok,-16:], interpolation = 'nearest', aspect = 0.1)
Out[9]:
In [10]:
TMPrimaryHeader = BitStruct('transfer_frame_version_number' / BitsInteger(2),
'spacecraft_id' / BitsInteger(10),
'virtual_channel_id' / BitsInteger(3),
'ocf_flag' / Flag,
'master_channel_frame_count' / BitsInteger(8),
'virtual_channel_frame_count' / BitsInteger(8),
'transfer_frame_secondary_header_flag' / Flag,
'synch_flag' / Flag,
'packet_order' / Flag,
'segment_length_id' / BitsInteger(2),
'first_header_pointer' / BitsInteger(11))
TMSecondaryHeaderID = BitStruct('transfer_frame_secondary_version_number' / BitsInteger(2),
'transfer_frame_secondary_header_length' / BitsInteger(6))
TMSecondaryHeader = Struct('id' / TMSecondaryHeaderID,
'data' / Bytes(this.id.transfer_frame_secondary_header_length + 1))
primary_headers = [TMPrimaryHeader.parse(bytes(f)) for f in frames[crc_ok]]
mcfc = np.array([p.master_channel_frame_count for p in primary_headers])
vcid = np.array([p.virtual_channel_id for p in primary_headers])
vcfc = np.array([p.virtual_channel_frame_count for p in primary_headers])
In [11]:
plt.figure(figsize = (8,4), facecolor = 'w')
plt.plot(mcfc)
plt.title('Master channel frame count')
plt.ylabel('Value')
plt.xlabel('Frame');
In [12]:
plt.figure(figsize = (8,4), facecolor = 'w')
plt.plot(np.diff(mcfc.astype('uint8')))
plt.title('Master channel frame count derivative')
plt.ylabel('Value')
plt.xlabel('Frame');
In [13]:
np.sum(np.diff(mcfc.astype('uint8')) - 1)
Out[13]:
In [14]:
vcids = np.unique(vcid)
for v in vcids:
print(f'Virtual channel ID {v}: {np.sum(vcid == v)} frames')
In [15]:
for v in vcids:
plt.figure(figsize = (8,4), facecolor = 'w')
style = '.-' if np.sum(vcid == v) < 100 else '-'
plt.plot(vcfc[vcid == v], style)
plt.title(f'Virtual channel {v} frame count')
plt.ylabel('Value')
plt.xlabel('Frame')
plt.figure(figsize = (8,4), facecolor = 'w')
plt.plot(np.diff(vcfc[vcid == v].astype('uint8')), style)
plt.title(f'Virtual channel {v} frame count derivative')
plt.ylabel('Value')
plt.xlabel('Frame')
In [16]:
for v in vcids:
plt.figure(figsize = (15,15), facecolor = 'w')
total = np.sum(vcid == v)
plt.imshow(frames[crc_ok][vcid == v], interpolation = 'nearest', aspect = 500/total)
plt.title(f'Virtual channel {v}')
The first header pointer in virtual channel 7 frames is always 2046, which means idle data only.
In [17]:
np.unique([p.first_header_pointer for p in primary_headers if p.virtual_channel_id == 7])
Out[17]:
This is a typical virtual channel 7 frame header:
In [18]:
[p for p in primary_headers if p.virtual_channel_id == 7][0]
Out[18]:
We now show the nature of the payload in these idle frames.
First we extract payload bits in each frame in virtual channel 7, leaving gaps (as zeros) for missing frames.
In [19]:
seq_unwrap = np.int32(np.round(np.unwrap(vcfc[vcid == 7]/256*2*np.pi)/(2*np.pi)*256))
vc7_bits = np.zeros((seq_unwrap[-1]-seq_unwrap[0]+1, (frame_size - 16)*8), dtype = 'float32')
for j, frame in enumerate(frames[crc_ok][vcid == v]):
vc7_bits[seq_unwrap[j]-seq_unwrap[0], :] = 2*np.unpackbits(frame[10:-6]).astype('float32')-1
Circular correlation analysis with FFT.
In [20]:
pn_len = 511
vc7_pn = vc7_bits.ravel()
vc7_pn = vc7_pn[:vc7_pn.size // pn_len * pn_len].reshape((-1,pn_len))
f = np.fft.fft(vc7_pn)
vc7_pn_corr = np.abs(np.fft.ifft(f * f[0,:].conjugate()))
vc7_pn_shift = np.argmax(vc7_pn_corr, axis = 1)
vc7_pn_corr = np.max(vc7_pn_corr, axis = 1)
In [21]:
plt.figure(figsize = (10,6), facecolor = 'w')
plt.plot(vc7_pn_shift)
plt.xlabel('Block')
plt.ylabel('Circular shift(bits)', color = 'C0')
ax2 = plt.gca().twinx()
ax2.plot(vc7_pn_corr, '.', color = 'C1')
plt.title(f'{pn_len} bit blocks circular correlation')
ax2.set_ylabel('Correlation magnitude (bits)', color = 'C1');
A more in detail look at the jumps in the correlation shift, cleaning the problems due to missing frames:
In [22]:
plt.plot(vc7_pn_shift[vc7_pn_corr>200])
Out[22]:
In [23]:
d_shift = np.diff(vc7_pn_shift[vc7_pn_corr>200])
d_shift[d_shift != 0]
Out[23]:
In [24]:
308 - 511
Out[24]:
In [25]:
idx_jumps = np.where(np.diff(vc7_pn_shift) % 511 == 308)
np.diff(idx_jumps)
Out[25]:
This motivates (and actually proves) the idea that the PN sequence is reset every 2250752 bits.
In [26]:
long_period = 2250752
We now build the same PN sequence, including the resets.
The PN sequence is generated by an LFSR with the polynomial x^9 + x^5 + 1.
In [27]:
lfsr = np.ones(9, dtype = 'uint8')
lfsr_out = np.empty(511, dtype = 'uint8')
for j in range(lfsr_out.size):
out = lfsr[8] ^ lfsr[4]
lfsr_out[j] = lfsr[8]
lfsr = np.roll(lfsr, 1)
lfsr[0] = out
This shows that the LFSR output (with an appropriate circular shift) coincides with the PN sequence.
In [28]:
lfsr_out_f = np.fft.fft(2*lfsr_out.astype('float32')-1)
lfsr_corr = np.fft.ifft(lfsr_out_f * f[0].conjugate()).real
np.argmax(lfsr_corr), np.max(lfsr_corr)
Out[28]:
We create our long LFSR output by repeating and truncating to 2250752 bits.
In [29]:
lfsr_out_long = np.tile(2*lfsr_out.astype('int')-1, int(np.ceil(long_period/511)))[:long_period]
We just need to repeat lfsr_out_long and trim it appropriately so that it lines up with the bits in virtual channel 7.
We know that we need to take trim away 294 bits due to the circular correlation offset, but we also need to find how many integer multiples of 511 bits we need to trim in addition. We derive this from the location where the first PN sequence reset happens in virtual channel 7.
In [30]:
int(np.ceil(long_period/511)) - idx_jumps[0][0] - 2
Out[30]:
So we trim 1806 full 511 bit periods.
In [31]:
lfsr_out_full = np.tile(lfsr_out_long, int(np.ceil(vc7_bits.size/lfsr_out_long.size)))[1806*511 + 294:]
Our LFSR resetting sequence only has -1's and 1's. The vc7_bits sequence has -1's and 1's, and 0's where we lost a frame. To check if they coincide (except for lost frames), we see that the product of both is never -1.
In [32]:
np.unique(lfsr_out_full)
Out[32]:
In [33]:
np.unique(vc7_bits.ravel())
Out[33]:
In [34]:
np.any(vc7_bits.ravel() * lfsr_out_full[:vc7_bits.size] == -1)
Out[34]:
This proves the structure of the idle data in virtual channel 7: A 511 bit PN sequence generated by the polynomial x^9 + x^5 + 1 with initial state 0x1ff and resetting every 2250752 bits.
It turns out that 2250752 bits is exactly 256 frames. Let us see if the sequence resets exactly when the virtual channel frame counter wraps to zero. It suffices to check if the frames with virtual channel frame counter equal to zero have a payload which starts by the PN sequence at shift zero.
In [35]:
wrap_frames = frames[crc_ok][(vcid == v) & (vcfc == 0)]
wrap_frames_start = np.unpackbits(wrap_frames[:,10:], axis = 1)[:,:511]
np.any(lfsr_out ^ wrap_frames_start)
Out[35]:
Indeed the sequence resets anytime that the virtual channel frame counter is zero.